Submitted to ApJ 

Preprint typeset using I4TgX style emulateapj v. 1 1/10/09 



THE DEBRIS DISK OF VEGA: A STEADY-STATE COLLISIONAL CASCADE, NATURALLY 

S. MULLER 

Astrophysikalisches Institut und Universitatssternwarte, Friedrich-Schiller-Universitat Jena, 
SchillergaBchen 2-3, 07745 Jena, Germany 



OS 

o 
o 

(N 

O 
CD 

Q 



PL, 

6 



> 
o 
o> 



(N 

o> 
o 



X 



T. LOHNE 

Astrophysikalisches Institut und Universitatssternwarte, Friedrich-Schiller-Universitat Jena, 
SchillergaBchen 2-3, 07745 Jena, Germany 

A. V. KRIVOV 

Astrophysikalisches Institut und Universitatssternwarte, Friedrich-Schiller-Universitat Jena, 
SchillergaBchen 2-3, 07745 Jena, Germany 

(Received 2009 September 7; Accepted 2009 November 30) 
Submitted to ApJ 

ABSTRACT 

It has been argued that the photometric data and images of the archetypical debris disk around Vega may 
be in contradiction with the standard, steady-state collisional scenario of the disk evolution. Here we perform 
physical modeling of the Vega disk "from the sources". We assume that dust is maintained by a "Kuiper belt" 
of parent planetesimals at ~ 100 AU and employ our collisional and radiative transfer codes to consistently 
model the size and radial distribution of the disk material and then thermal emission of dust. In doing so, we 
vary a broad set of parameters, including the stellar properties, the exact location, extension, and dynamical 
excitation of the planetesimal belt, chemical composition of solids, and the collisional prescription. We are 
able to reproduce the spectral energy distribution in the entire wavelength range from the near-infrared to 
millimeter, as well as the mid-IR and sub-millimeter radial brightness profiles of the Vega disk. Thus our 
results suggest that the Vega disk observations are compatible with a steady-state collisional dust production, 
and we put important constraints on the disk parameters and physical processes that sustain it. The total disk 
mass in < 100 km-sized bodies is estimated to be ~ 10 Earth masses. Provided that collisional cascade has 
been operating over much of the Vega age of ~ 350 Myr, the disk must have lost a few Earth masses of solids 
during that time. We also demonstrate that using an intermediate luminosity of the star between the pole and 
the equator, as derived from its fast rotation, is required to reproduce the debris disk observations. Finally, we 
show that including cratering collisions into the model is mandatory. 

Subject headings: planetary systems: formation - circumstellar matter - radiation mechanisms: thermal - stars: 
individuals: Vega 



1. INTRODUCTION 

Debris disks - disks of planetesimals and dust around main- 
sequence stars - are known to encircle a sizeable fraction of 
main-sequence (MS) stars. The observed dust is believed to 
derive from collisions between unobserved planetesimals that 
were neither used to make up planets nor ejected from the 
system by th e time when the nebular gas was dispersed (see 
IWyattll2008[ for a review). 

A standard scenario of a debris disk evolution (e.g. 
Thebault et al.1 120031: iKrivov et all l2006t IStrubbe & Chiang! 
2006t iThebault & Augereaul 120071: lLohneetal.1 120081: 



Krivov et al. 2008; Thebault & Wu 2008) can be summarized 



as follows. There is a relatively narrow belt of planetesimals 
("birth ring") in orbits with moderate eccentricities and 
inclinations, exemplified by the classical Kuiper belt in the 
solar system. The orbiting planetesimals in the birth ring 
undergo collisional cascade that grinds the solids down 
to dust. At smallest dust sizes, stellar radiation pressure 
effectively reduces the mass of the central star and quickly 
(on the dynamical timescale) sends the grains into more 
eccentric orbits, with their pericenters still residing within 
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the birth ring while the apocenters are located outside the 
ring. As a result, the dust disk extends outward from the 
planetesimal belt. The smaller the grains, the more extended 
their "partial" disk. The tiniest dust grains, for which the 
radiation pressure effectively reduces the physical mass by 
half, are blown out of the system in hyperbolic orbits. The 
radiation pressure blowout of the smallest collisional debris 
represents the main mass loss channel in such a disk. A 
steady state implies a balance between the production of dust 
by the collisional cascade and its losses by radiation pressure 
blowout. 

Such disks are usually referred to as collision-dominated, 
as opposed to s ystems that — at du st sizes — are transport- 
dominated (e.g. IKrivov et al.l 120001) . In the latter case, ra- 
dial transport of dust material by various drag forces occurs 
on shorter timescales then collisions. Then, additional re- 
moval mechanisms may play a significant role. For example, 
Poynting-Robertson (P-R) drag can bring grains close to the 
star where they would sublimate or deliver them into the plan- 
etary region where they would be scattered by planets. To be 
transport-dominated, the system should either have an opti- 
cal depth below current detection limits ( Wyatt 2005a) or be 
subject to transport mechanisms other than P-R drag, such as 
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strong stellar winds around late -type stars (Strubbe & Chiang] 
l2006tlAugereau & BeusT2 006). Most of debris disks detected 
so far are thought to be collision-dominated. Accordingly, 
transport-dominated systems are not considered here. 

The size segregation in collision-dominated disks de- 
scribed above means that at different wavelengths the 
same disk would look d ifferently (see, e.g., Fig. 17 in 
Theba ult & Augereaull2007l) . Measurements at longer wave- 
lengths (sub-mm) probe larger grains, because they are cooler, 
and thus trace the parent ring. Such observations may also re- 
veal clumps, if for instance there is a planet just interior to the 
inner rim of the parent ring, and parent planetesimals and their 
debris are trapped in external reso nances, similar to Plutinos 
in 3 : 2 resonance with Neptune dWvatlll2~006t iRrivov et al.l 
2007). At shorter wavelengths (far-IR, mid-IR), smaller 
(warmer) grains are probed. Thus the same disk appears much 
larger. It may appear featureless, even if t he parent rin g is 
clumpy, because strong radiation pressure jW yatt 200g) and 
non-negligible relative velocities dKrivov et al.l 20071) would 
liberate such particles from resonant clumps. As a result, they 
would form an extended disk, as described above, regardless 
of whether their parent bodies are resonant or not. 

The observed statistics of infrared (IR) luminosities of de- 
bris disks around AFGK stars, including their dependence on 
the stellar age and spectral class, is largely consistent with 
the steady-state coll i sional evolution sce nario presented above 
dWvatt et al.M2007bh Eohne et all 120081) . Furthermore, a de- 
tailed spectral energy distribution (SED) from far-IR to mm 
wavelengths, which was measured for a handful stars, can be 
repro duced with the mod els of a steady-state collisional cas- 
cade (Krivov et al. 2008). There are some exceptions, how- 
ever (see. e.g.. lWvattll2008L and references therein). For in- 
stance, some A-type stars have fractional luminosities that 
are too high compared to what is expected from collisionally 
evolving Kuiper belt analogs. A few percent of FGK stars 
show excess emission shortward of 24 /mi that must come 
from "asteroidal" region inside ~ 10 AU, which in some cases 
is also too high to b e compatible with the standard scenario 
dWvatt et al.ll2007al) . Besides, incidences and properties of 
debris disks around M-type stars remain uncertain: the obser- 
vations are scarce, and the physics of such disks may be dif- 
feren t due to low stellar luminosity and strong stellar winds 
(e.g. JStrubbe & Chiang|l2006t lAugereau & B eust 2006|. De- 
spite these caveats, the majority of the debris disks discovered 
so far appears compatible with the steady-state collisional sce- 
nario. However, the vast majority of debris disks are as yet un- 
resolved, and the set of observables is typically limited to the 
fractional luminosity and a few photometric points in the IR. 
Many more constraints on the disk properties could be posed 
by resolved images, if these available, especially at more than 
one wavelength. There are currently a dozen of disks with 
such datasets. Accordingly, in this paper we undertake prob- 
ably one of the first attempts to investigate, whether all avail- 
able data on one particular debris disk star are compatible with 
a steady-state collisional scenario of dust production and evo- 
lution. We choose Vega. 

Vega was the first MS star around which an IR excess 
over the ph otospheric flux, indi cative of debris dust, was 
discovered (lAumann et alJ 119841) . It is therefore treated 
as an archetypical debris disk star; MS stars with IR 
excesses are often named "Vega-type stars". Fluxes at 
wavelengths fro m mid-IR to millimeter have been mea- 
sured with I RAS (lAumann et al.|1984UWalker & Wolstencrofj 
119881) . KAO dHarper et alJll984HHarvev et al.ll 19841) . and ISO 



dHeinrichsen et al.|[l9 98). As a result, its SED is known rela- 
tively well. The Vega disk has been resolved in sub-mm and 
mm w ith JC MT dZuckerman & Becklinlll993t iHolland et all 
1 19981). OVRO dKoerner et al.l l2001h. IRAM dChini et al.l l 19901: 
Iwilner et all 120021) . and CSO dMarsh et all 120061 These 
observations reveal a clumpy ring of large dust grains be- 
twee n abou t 80 an d 12 AU, suggesting a Kuiper belt ana- 
log. IWvattl d2003l) and lReche eTaLl (12008b naturally explain 
the ring structure with a resonant trapping of dust parent bod- 
ies by a presumed Neptune- to Saturn-mass planet during their 
o utward m i gratio n in the past. 

ISu et all (120051) resolved the Vega system by Spitzer/MIPS 
at 24, 70, 160jum. They found a featureless, huge disk ex- 
tending up to ~ 800 AU. Although it came as a surprise at the 
time when this discovery was made, it is no longer astonish- 
ing now. As noted above, this is exactly what is expected: a 
Kuiper belt-sized, clumpy ring of large dust grains seen in the 
sub-mm and a much more extended disk of small grains, pro- 
ducing a smooth brightness distribution evi dent in the mid- 
to far-IR. However, by fitting these data, ISu et all (12005b 
deduced a mysterious overabundance of blowout grains of 
~ 1 /mi in radius. Under the assumption of a steady-state col- 
lisional disk evolution over the Vega age, this would imply 
that the disk must have lost ~ 3 Jupiter masses of material, 
which appeared unlikely. Accordingly, they suggested a re- 
cent major collisional event as a possible explanation. More 
exotic alternative scenarios proposed to explain such a large 
fraction of blowout g rains include a close stellar encounter 
dMakarov et alJl2005l) and a dynamical instability event sim- 
ilar to wha t caused the Late H eavy Bo mbardment in the so- 
lar sys tem dWvatt et al.l2007aT) . However, Kenyon & Bromley 
(2008), who modeled the Vega debris disk as an aftermath of 
icy planet formation with their hybrid multi-annulus coagula- 
tion code, find their model to be capable of reproducing the 
Spitzer fluxes, questioning the need in alternative scenarios. 

We note that the excessive a mount of grains in blowout or- 
bits inferred bv lSu et ail (2005) uncovers another problem. A 
steady-state collisional evolution implies a certain size distri- 
bution of dust. Typically, the amount of blowout grains in- 
stantaneously present in the steady-state system is much less 
than the amount of slightly larger grains in loosely bound or- 
bits around the star. This is because the dust production of 
the grains of adjacent sizes is comparable, but the lifetime of 
bound grains (due to collisions) is much longer than the life- 
time of blowout grains (disk-crossing timescale). A jump in 
the size distribution around the blowout size is a robust pre- 
diction of all collisional cascade models, even those that do 
not assume a steady state (e.g. those that describe short-lived 
consequences of a major brea k-up event). T hus the amount 
of blowout grains reported by Su et al. (2005) would necessi- 
tate an unrealistically huge amount of larger grains in bound 
orbits. And this conclusion would hold not only in a steady- 
state scenario, but also all alternative scenarios listed above 
would face the same problem. 

Beside the outer disk, the inner part of the system reveals 
another peculiarity. Pioneering inter ferometry observ ations 
with CHARA/FLUOR in the near-IR dAbsil et alj|2006h have 
led to the discovery of a dust cloud just exterior of the subli- 
mation zone, well inside 1 AU ("exozodi"). 

Just like the d ust in the system, the c entra l star turned out to 
be unu sual, too. iPeterson et alJ d2006h and lAufdenberg et al] 
(2006) found Vega to be a rapid rotator, which makes stellar 
parameters functions of the stellar latitude. Table Q] summa- 
rizes the stellar parameters relevant for this study. It remains 
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TABLE 1 
Stellar parameters. 





Equator Pole 


Note 


R[R e ] 


2.873 ±0.026 2.306 ±0.031 


1 


T eS [K] 


7 900+50° 10 150 ±100 


2 


U [is] 


28 + j? 57 ± 3 


3 


log(g[cms~ 2 ]) 


4.074 ±0.012 3.589 ±0.056 


1 


M, [M ] 


2.3 ± 0.2 


2 


Age [Myr] 


350 





Notes : (1) From [Peterson et al. (2006), (2) From Aufdenbe rg et alj 
(2006), (3) Luminosity at the equator and at the poles derived from 
the equatorial and polar values of the stellar radius R and effective 
temperature T^r an d from the average stellar luminosity of 37 ± 3 L Q 
(Aufdenberg et al. 2006) through the Stefan-Boltzmann relation. 

unclear whether unusual properties of the disk are somehow 
related to those of the star. 

In this paper, we re-address the question of whether the 
available data are compatible with a steady-state collisional 
scenario of dust production and evolution. Instead of simply 
seeking dust distributions - e.g. in the form of power laws 
- that would provide the best fit to the observables , we em- 
ploy an approach described by iRrivov et alj (120081) . In this 
approach, we assume a planetesimal belt with certain prop- 
erties, evolve it with a collisional code to generate the dust 
portion of the debris disk, calculate the emission of that dust, 
and compare it to the observed emission. The best fit can be 
achieved by varying the parameters that describe the planetes- 
imal belt (e.g. its location and mass), rather than parameters 
of the dust distribution as is commonly done. 

In Sect. [2] we describe the data and their reduction and in 
Sect. [3] dynamical and thermal emission models used in this 
paper. Section [4] presents our reference ("first-guess") disk 
model and compares it with the available observational data. 
In Sect. [5] we analyze the influence of various model parame- 
ters on the observables. Section [6] contains a discussion and 
Sect.|7]lists our conclusions. 

2. OBSERVATIONAL DATA 

From the Spitzer archive, we extracted MIPS images of the 
Vega system at 24, 70, and 160 /mi, using the Leopard soft- 
warq']. After applying standard basic corrections, the data 
were further processed to remove remaining constant back- 
grounds. Column averages were used to eliminate the so- 
called jail-bar artifacts caused by saturated sources on 24 pm 
images. 

For each image, the center of the source was then deter- 
mined by minimizing first order moments. For the 160//m 
images, the pointing information was used for that purpose. 
After that, radial profiles were derived in an iterative process 
with two steps: (i) sampling the interpolated image at dis- 
cretely binned distances around the source center with a sub- 
pixel resolution and (ii) integrating and subtracting the result- 
ing profile pixel-wise from the original. Five such iterations 
were performed for each image. Finally, all profiles for each 
wavelength were combined, using their median. 

At 24 pm, the data reduction is complicated by the fact that 
the central part of the images (< 4" or < 30 AU) is satu- 
rated. Still, after subtracting the reference Point Spread Func- 
tion (PSF), the outer part of the brightness profile can be found 

1 http://ssc.spitzer.caltech.edu/mips/ 
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FIG. 1. — Radial profiles of the surface brightness for Vega as extracted 
from the Spitzer MIPS images for 24 (large squares), 70 (large circles), 
and 160//m (large triangles) from top to bottom. The standard sources 
HD 217382, Sirius, and (40) Harmonia used for PSF subtraction for the dif- 
ferent wavelengths are shown as small symbols of the same shapes. 

with sufficient accuracy. Figure [T](top) shows the median pro- 
file extracted for Vega compared to that of the reference star 
HD 217382, a K4III giant at a distance of 120 AU, sca led 
up from its intrinsic 2.15 Jy to the 7.4 Jy dSu et al.1120051) of 
Vega's photosphere. The resulting photometric excess, ob- 
tained by integrating the subtracted profile (the Vega profile 
minus the scaled PSF) from 4" outward, is 0.94 + 0.28 Jy. 
Here the formal error is determined by the standard deviations 
of Vega and the PSF plus an assumed systematic error of 5 % 
in the photometry of Vega's photosphere. However, the flux 
is probably more uncertain because so is the contribution of 
the inner part of the disk, between 4" and 10" (30-80 AU). 
Thus, for a more accurate comparison with our models below, 
we can use the "certain" part of the observed flux by inte- 
grating the brightness profile from 10" outward. This nearly 
halves the total 24 pm flux, giving 0.53 Jy. This "partial" flux 
will later be compared with the flux predicted by our models 
exactly in the same range of distances, from 10" (80 AU) out- 
ward. Coincidentally, it is this range where the emitting dust 
is only present in most of our models, because the inner edge 
of the birth ring is located at 80 AU. We emphasize, however, 
that the true total 24 pm excess flux may be higher. It proba- 
bly lies in the range 0.5-0.9 Jy, with no obvious possibility to 
narrow it because of the saturation problem described above. 

The reduction of the 70 pm images was more straightfor- 
ward, resulting in the radial profiles shown in Fig.Q~](middle) 
and an excess of 8.6 + 1.2 Jy. Sirius (HD 48915) was used as 
reference. 

Since stray light strongly contaminates images of hot, 
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TABLE 2 
Photometric data. 



A [jiml F dlsk [mJy] Instr. 

25.0 (2.979 ± 0.936) x 10-' IRAS 

25.0 (4.019 ± 2.080) x 10 3 ISO 

47.0 (6.414 ± 1.640) x 10 3 KAO 

60.0 (7.918 ± 0.901) x 1() 3 IRAS 

60.0 (9.318 ± 2.082) x 10 3 ISO 

80.0 (9.091 ±2.910) x 10 3 ISO 

95.0 (6.829 ± 1.825) x 10 3 KAO 

100.0 (7.109 ± 0.754) x 10 3 IRAS 

100.0 (5.969 ± 1.920) x 10 3 ISO 

120.0 (4.890 ± 1.551) x 10 3 ISO 

170.0 (2.437 ± 0.771) x 10 3 ISO 

193.0 (8.932 ± 5.000) x 1() 2 KAO 

200.0 (1.321 ±0.408) x 10 3 ISO 

350.0 (4.691 ± 1.500) x 10 2 CSO 

450.0 (1.301 ±0.450) x 10 2 CSO 

800.0 (1.626 ± 0.500) x 10' IRAM 

800.0 (2.426 ± 1.500) x 10' JCMT 

850.0 (4.086 ± 0.540) x 10' JCMT 

870.0 (1.721 ±0.900) x 10' IRAM 

1300.0 (2.390 ± 1.500) x 10° IRAM 

1300.0 (1.189 ± 0.190) x 10' OVRO 

1300.0 (1.140 ± 0.170) x 10 1 IRAM 



Ref. 



Walker & Wolstencroft (1988) 
Heinrichsen et al. (1998) 
Harvey et al. (1984) 
Walker & Wolstencroft (1988) 
Heinrichsen et al. (1998) 
Heinrichsen et al. (1998) 
Harvey et al. (1984) 
Walker & Wolstencroft (1988) 
Heinrichsen et al. (1998) 
Heinrichsen et al. (1998) 
Heinrichsen et al. (1998) 
Harper et al. (1984) 
Heinrichsen et al. (1998) 
Marsh et al. (2006) 
Marsh et al. (2006) 
Chini etal. (1990) 
Zuckerman & Becklin ( 
Holland etal. (1998) 
Chini et al. (1990) 
Chini etal. (1990) 
Koerneretal. (2001) 
Wilner et al. (2002) 



i.e. stellar, sources a t 160/vm, we followed the strategy of 
Stape lfeldt et alJ (l2004f) and used an additional reference star 
(HD 197989) to remove this artifact. However, visual inspec- 
tion of the Vega images already suggested that the removal 
procedure would not be perfect because of some dissimilar- 
ities in the artifact shapes. Therefore, only the half-image 
facing away from the a rtifact center w as used for profile ex- 
traction, as was done bv lSu et alJ ( 120 05). The actual reference 
used for PSF subtraction was the asteroid (40) Harmonia, al- 
though its influence on the excess of 1.9*9 ! Jy is negligible. 
See Fig. Q] (bottom) for the resulting profiles. 

The photometry points across the wavelength range from 
mid-IR to millimeter taken from the literature, are listed in 
Table [2] These points are plotted as symbols with error bars 
in bottom left panels of Figs. |5hT3] 

3. METHODS 
3.1. Dynamical and Collisional Evolution 

Our technique to follow the size and radial distribution 
of material in rotationally-symm etric debris disks is de- 
scribed in detail in previous pa pers (Kri vov et al.| [2000, 2005, 
l2006l l200l lLohne etal.1120081) . Our numerical code, ACE 
(Analysis of Collisional Evolution), solves the Boltzmann- 
Smoluchowski kinetic equation over a grid of masses, perias- 
tron distances, and orbital eccentricities of solids. It includes 
the effects of stellar gravity, direct radiation pressure, drag 
forces, as well as disruptive and erosive collisions. 

The results of ACE simulations depend sensitively on the 
adopted collisional prescription. Possible collisional out- 
comes range from a perfect sticking or fragmentation with 
subsequent reaccumulation of particles (when the impact en- 
ergy is low) to a crater formation or even a complete disrup- 
tion (when it is high). An important quantity in the colli- 
sional prescription is the critical specific energy for disrup- 
tion and dispersal, Q^. It is defined as the impact energy per 
unit target mass that results in the largest remnant containing 
a half of the original target mass. For small objects, is 
determined solely by the material strength, while for objects 



larger than ~ 100 m, the gravitational binding energy domi- 
nates. As a result, is commonly described by the sum of 
two power laws (see, e .g.. iDavis et al][T9 85; Holsappl ell 19941: 
Paolicchi et alJ [19961: iDu rda & Dermo tj Il997t iDurdaetal] 
1998tlBenz & Asphaug||1999HKenvon & Bro mlev 20043): 



Qo = Qd,s(s)(-^) +Qd, b (s)(j 



-) 



3b, 



km/ 



(1) 



where the subscripts s and g denote the strength and gravity 
regime, respectively. 

The thermal emission in the mid- and far-IR is dominated 
by small dust with sizes slightly above the so-called blowout 
limit, below which particles are repelled from the system by 
stellar radiation pressure. An equilibrium size distribution set 
by the collisional cascade, especially near that blowout limit, 
depends on the total mass and the size distribution of frag- 
ments of each individual collision. Depending on the ratio of 
the impact energy to the critical energy, the avera ge size of the 
fragm ents changes. This has been measured (e.g. Taka gi et al.1 
1984) and numerically modeled (e.g. lDurda et al.ll2007l) . and 
useful analytic prescriptions for use in co llisional codes are 
available (e.g. Thebault & Augereau 200?]). Here we apply a 
simplified approach with only one parameter, the slope rj of a 
single power law AN oc m~ n dm. Larger values of r\ put more 
weight on small fragments, while smaller values of r\ "prefer" 
larger fragments. Note that using a single power law down to 
infinitely small grains is only meaningful for r\ < 2. 

Each ACE run outputs, among other quantities, the size and 
radial distribution of disk solids over a broad size range from 
sub-micrometers to hundreds of kilometers at different time 
steps, and the code is fast enough to evolve the distribution 
over gigayears. As shown in Fig. [21 typical timescales for P- 
R drag in the Vega disk are much longer than collisional life- 
times, except for a very narrow size range close to the blowout 
limit, where both become comparable. Thus we switch off 
the P-R effe ct in the ACE runs, but make additional checks 
in Sect. 16.81 Gas drag can safely be neglected, because the 
Vega system with its 350 Myr age could not have retained any 
primordial gas, while the density of secondary gas cannot be 
high enough to influence the dust dynamics. As long as the 
drag forces are discarded, the number of parameters to vary, 
and thus the number of required A CE runs, can be re duced 
by applying scaling laws derived in lLohne et"aTI (120081) and 
iKrivov et al.1 J2008). Specifically, it is not needed to perform 
separate ACE runs for the initial planetesimal disk with dif- 
ferent initial masses. 

3.2. Thermal Emission of Dust 

Considering a rotationally symmetric disk of spherical par- 
ticles, we denote by N(r, s) the disk's surface number density 
of grains of radius s at a distance r from the star. The specific 
thermal emission flux from the entire disk Ff% sV measured 
by an observer at a distance D from the star at a given wave- 
length A, can be calculated as (Kri vov et alJ 12008) 



r /Ldisk' 



2^ 
D 2 



I 



dr(T„) 



dr P 



I 



ds s x 



xN(r,s)Qf s (s)B A (TX 



(2) 



where B^Tg) is the Planck function, Q a A (s) is the absorption 
efficiency, A is the wavelength, and the grain temperatures 
T g (s) are deduced from the standard assumption of thermal 
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FIG. 2. — Timescales as functions of the grain size: collisional time (solid 
line), blowout time (long-dashed), and P-R time (short-dashed). The colli- 
sional time is an average over the grain orbits with all possible pericentric 
distances q and eccentricities e, weighted with the amounts of those parti- 
cles in the disk. The P-R time is the time it takes for a grain to drift across 
the parent ring or, more exactly, the time interval over which a grain's q re- 
duces from 120 AU to 80 AU. It was computed by simul taneously solvin g the 
orbit-averaged equations that describe q(i) and e(t) (Burns et al. 1979). All 
timescales are for the parameters of the reference model (see Sect.|4). 

equilibrium with the stellar radiation field. Similarly, the ra- 
dial profile of the surface brightness S ^(r) is given by 

SA(r) ^^j ds s2 N(r, s)Qf s (s)B A (T g (r, s)). (3) 

To enable comparison with observations, equation (f3]l has to 
be convolved with the normalized PSF of the instrument: 

5»=5^(r)®PSF A (r). (4) 

A numerical solution of Eqs. <J2J and ^ has been imple- 
mented in two routines, SEDUCE {SED Utility for Circum- 
stellar Environment) and SUBITO (SUrface Brightness Inves- 
tigation TOol), respectively. The stellar flux needed to com- 
pute the dust temperatu res is extracted from NextGen models 
dHauschildt et al.ll 19991) . Both codes have a direct interface to 
ACE, so that ACE output can be used as input to SEDUCE 
and SUBITO. To arrive at the correct total dust mass and ab- 
solute values of fluxes, we rescale the ACE-SEDUCE and 
ACE-SUBITO simulation results with the aid of t he appropri- 
ate sca ling laws, as described in Appendix A of iRrivov et al.l 
(l2008h . 

4. THE REFERENCE MODEL OF THE VEGA DISK 

4. 1 . Choice of Model Parameters 

In the reference model, we use the stellar luminosity at the 
equator of L+ = 28 L Q (Table [TJ, and we assume that the 
collisional cascade has been operating over the entire stellar 
age, 350 Myr. 

According to iDentetai] d2000b ISu et alj d2005l) . 
iMarsh et all (120061) and IWyattl (120061) we adopt — as a 
"first guess" — an initial ring of parent bodies with semima- 
jor axes ranging from 80 to 120 AU and an initially constant 
surface density in this range. The clumpy shape of the 
sub-mm ring, usually interpreted through resonant capture of 
planetesimals by an unseen planet interior to the ring, implies 
that the eccentricities of the plan etesimals are not very low 
(Wyatt 2001 iReche et all 120081) . On the other hand, the 
relatively narrow ring observed at wavelengths longer than 
350 j-im indicates that eccentricities cannot be too high. As a 
reasonable compromise and for the sake of simplicity, for the 
reference model we adopt a uniform distribution of eccentric- 
ities from 0.0 to a moderate value of 0.2. Maximum orbital 



inclinations (or a semi-opening angle) were then set to 0.1 rad 
in accord with the energy equipartition assumption. Thus, 
the initial planetesimal disk resides between 64 and 144 AU 
from the star. This is still in agreement with the observed 
80 to 120 AU as most of the material is concentrated in the 
central part of the initial ring. Note that these assumptions 
describe the initial disk extension. The subsequent collisional 
and dynamical evolution of the parent belt slightly changes 
the distributions of planetesimals. 

All particles were assu med to be composed o f astronomi- 
cal silicate (p = 3.3 gcm~ 3 . lLaor & Draindll993l) . Mie theory 
was used to calculate radiation pressure and absorption effi- 
ciencies (Fig.0. 

The disk was modeled with ACE with grains ranging from 
0.05 fim to 67 km in radius and the mass ratio in the ad- 
jacent bins of 4. The pericenter distance grid covered 50 
logarithmically-spaced values from 20 AU to 800 AU. The ec- 
centricity grid contained 100 linearly-spaced values between 
-5.0 and 5.0 (eccentricities are negative in the case of small- 
est grains with B > 1, whose orbits are anomalous hyper- 
bolas, open outward from the star). The distance grid used 
by ACE to output distance-dependent quantities such as the 
size distribution was 10 AU through 600 AU at 10 AU in- 
crements. In the collisional prescription, we set (2d,s(s) = 
gD, g O) = 5 x 10 6 ergg- 1 , 3b g = 1.38, and 3b s = 0.37 and 
take the size distribution of fragments to be a power law with 
index rj = 1.833. Both disruptive and cratering collisions are 
switched on. 

A set of parameters used for the reference model is given in 
the first line of Tab. [3] In the table, we only list those param- 
eters that we later vary with respect to the reference model. 

4.2. Size, Radial, and Temperature Distributions 

Dust distributions in our reference model are presented in 
Fig. E] The right panel shows the grain size distributions 
within and outside the birth ring. The radiation pressure 
blowout effect causes a steep drop between 3 and 5 //m, which 
corresponds to B a 0.5 (Fig. [3] right). As previous studies 
have shown, the blowout drop in the size distribution results 
in a more or less pronoun ced wavy pattern in the distribu- 
tion o f larger particles (e.g. lThebault et al.ll2003l:lKrivov et al.l 
2006),with the "wavelength" and amplitude of the pattern de- 
pending on material strength and i mpact velocities. Ho wever, 
compared to previous studies (e.g. IRrivov et al. 2006), ACE 
now uses a grid of pericentric distances instead of semima- 
jor axes. This reduces the effects of discretization on the ef- 
fective relative velocities, which is strongest for particles on 
highly eccentric orbits with pericenters close to the birth ring. 
Therefore, the waviness is washed out, especially in the refer- 
ence run. 

The radial distribution of different-sized grains is shown in 
the bottom panel of Fig. |4]with thin lines. As explained above, 
most of the material is located between 80 and 120 Alfl The 
largest particles are confined to this region as they are nearly 
unaffected by radiation pressure. The smaller the particles, 
the wider they are spread over the disk. In addition, the bot- 
tom panel in Fig. |4]plots the total normal geometrical optical 

2 Note that the ring stalls at about 73 AU and not at the above-mentioned 
64 AU. This is due to the eccentricity binning. Individual bins are 0. 1 wide 
and centered at 0.05, 0.15, and so on. Thus, the largest effective eccentricity 
for e max = 0.2 is e = 0.15. The corresponding minimum pericentric distance 
for a = 80 AU is q min = 0.85 X 80 AU = 68AU. The next larger point in the 
pericenter grid is then centered at 73 AU. 
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FIG. 3. — Left: Absorption efficiency for 1 (thin lines) and 4.9/im (thick lines) particles consisting of astronomical silicate (solid lines) and of an astronomic 
silicate matrix with 10 % iron (dashed lines) and water ice (dotted lines) inclusions. Right: The ft ratio of the same grains as in the left panel (the same linestyles), 
assuming L = 28 L Q . Additionally, thin and thick dash-dotted lines show the /? ratio for pure astrosilicate grains, but higher stellar luminosities of L = 37 Lq and 
57 L Q respectively. 



TABLE 3 

Sets of model parameters used in the simulations. The parameters of the reference model are listed in the first line. For all 
other models, only parameters that are different from those of the reference model are given. 



Model L, [L ] T [Myr] a macr [AV] Q ou ter [AU] e max W composition collisions go.., [ergg '] b s rj 



ref. 28 350 80 120 0.2 0.1 no incl. w/cratering 5.0 xlO 6 0.37 1.833 

al — 35 — ____ _ ___ 

a2 — 3.5 — ____ _ 

b 50 ____ _ ___ 

b2 — — 100 ____ _ ___ 

b3 — — — 100 — — — — 

M — 150 — — — — ___ 

cl 37 — — ____ _ ___ 

c2 57 — — ____ _ ___ 

c3 — — — — — — ice incl. — — — — 

c4 — — — — — — iron incl. — — — — 

dl — — 71.1 130.9 0.1 0.05 — — — — — 

d2 — — 91.4 100.8 0.3 0.15 — — — — — 

el — — — — — — — w/o cratering — — — 

e2 — — — — — — — — 2.5 x 10 6 0.2 — 

e3 — — — — — — — — 6.9 x 10 6 0.45 — 

fl — — — ____ _ — — 1.6 

£2 — — — ____ _ __ 1.95 

"fit 45 = 62 120 01 O05 = = = = 1.95 



depth r (in arbitrary scale). The optical depth beyond about 
120 AU is dominated by particles just above the blowout size, 
~ 5 /an, which are in barely bound orbits. It is only in the 
region of the birth ring where particles with s > 10 /mi make 
a significant contribution to the optical depth. 

To judge about the dust temperatures and the thermal emis- 
sion of the disk in the reference model, the left top panel in 
Fig. |4] depicts the dust temperature as a function of stellar 
distance and grain size. Both distance and size axes share 
those of the two other panels. Hence, extrapolating the max- 
ima of the size and radial distributions into the temperature 
plot yields a "typical" temperature, i.e. the temperature of 
grains with r-dominating sizes at the distance where r peaks. 
For our reference disk these are grains with s w 5... 10 /mi at 
r ~ 80... 120 AU and their temperature is about 60 K. 

4.3. Dust Mass and Disk Mass 



We then used SEDUCE to obtain the SED or the refer- 
ence disk and to fit it vertically to the available observations 
starting at 25 /mi. Shorter wavelengths were neglected in the 
fitting process as the uncertainty of the photospheric sub- 
traction there is too high. Furthermore, we did not use the 
Spitzer 24 //m data point for fitting because of the uncertain- 
ties in converting the images into photometry, as discussed in 
Sect. O. In the SED calculation, we adopted the stellar pa- 
rameters at the equator to obtain the photosphere seen by the 
dust disk, but took the polar values to calculate the observed 
photosphere (Tab. Q]). 

The aforementioned fitting itself is not as straightforward 
as it may seem. As we wish the modeled absolute thermal 
emission flux to match the actually observed one, we need 
to change the amount of dust, i.e. the dust mass. In our 
approach, however, only the parameters of the parent plan- 
etesimals can be changed, not those of dust they produce. It 



The Debris Disk of Vega 



7 



200 



10- 



10" 



100 E 



50 NT 



10"'° 
10 17 

'§ 10 19 

"l 10-20 

< 10 21 
10 22 
10' 23 



24 urn 




'II 

^ i I 


: 


: 












* i/ 

- Size distribution at: ■' / 


8 


OAU > / 


100 AU ; / 


120 AU : / 


200 AU / 



total optical depth 



Radial distribution for: 

0.9 urn 

5.1 um --■ 
11 pjn -- 
96 urn — 



10" 10 _ " 10"" 10" 
A [cm 2 cm" 3 ] 



10 



100 



r[AU] 



10" 



E 

n 1 I 



10 3 



FIG. 4. — Left top: Disk's temperature profile for the assumed astrosilicate grains and stellar parameters as derived for Vega's equator. Two solid lines 
separate the regions of dominant emission at the three Spitzer/MIPS wavelengths: in the left part 24 /jm emission is more efficient than emission at the other two 
wavelengths, in the central part 70 //m emission dominates, and in the right part 1 60pm emission is the strongest. Left bottom: Radial distribution in the reference 
model for 0.9, 5.1, 11 and 96/jm grains (thin lines) and the resulting total optical depth in arbitrary units (thick solid line). Right: Grain size distribution in the 
reference disk at 80, 100, 120, and 200 AU. 



means that we have to go one step back and modify the ini- 
tial disk mass. However, it is not sufficient to change the ini- 
tial disk mass by the ratio of the observed and the modeled 
fluxes. The reason is that a change in the initial mass also al- 
ters the rate of the collisional evolution, whereas we need the 
"right" flux at a fixed time instant, namely — in the reference 
model — at 350 Myr. Therefore, to find the mass modifica- 
tio n factor we apply s caling rules, as explained in Appendix A 
of dKrivov et alj|2008l) . 

The dust, disk, and initial disk masses in the reference 
model are given in the first line of Tab. [4] The dust mass 
of 7 x 10~ 3 M m is by about a factor of two higher than what 
was derived bv lSu et all Q2006). The actual agreement is even 
better, because o ur upper cutoff size of 1000 /mi is larger than 
that of ISu et aTl The total mass of the reference disk is about 
16 M m , which is 85 % of its initial mass 350 Myr ago when 
the collisional cascade started to operate. 

4.4. Spectral Energy Distribution 

For an easier comparison between the modeled SEDs and 
the photometric observations in the different spectral regions, 
throughout the paper we use the excess ratio. The latter is 
defined as the ratio of the dust emission to the stellar photo- 
spheric emission or equivalently, as the ratio of the total flux 
(star + dust) to the stellar flux minus unity. The SED of the 
reference disk in terms of the excess ratio is presented in Fig. 
[5] (left bottom) with a solid line. Given that our reference 
model is a first-guess one, the agreement with the observa- 
tions is quite satisfactory. 

At 24 //m the model yields 0.43 Jy, which is at 2cr under 



TABLE 4 

Derived dust masses, disk masses, and initial disk masses for 
all models 



runs 


M dust [10- j M e ] 


M disk [M«] 


M u [M ffl ] 


ref. 


6.63 


16.3 


18.9 


al 


5.62 


4.21 


4.35 


a2 


3.67 


3.05 


3.06 


bl 


5.41 


19.3 


22.7 


b2 


7.34 


14.1 


16.2 


b3 


6.22 


17.1 


20.7 


b4 


7.38 


14.2 


15.7 


cl 


6.96 


14.8 


17.1 


c2 


6.59 


11.9 


13.6 


c3 


7.08 


18.5 


21.5 


c4 


7.28 


16.5 


19.1 


dl 


10.9 


41.5 


47.5 


d2 


4.74 


10.5 


12.6 


cl 


4.37 


3.60 


3.81 


c2 


8.64 


51.0 


62.7 


c3 


4.77 


5.15 


5.70 


fl 


7.09 


9.31 


10.5 


f2 


5.92 


32.2 


39.4 


fit 


5.09 


46.7 


55.5 
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FIG. 5. — Top left: Grain size dust distributions of the reference model (thick solid line) and of the same model but at earlier times of 3.5 (dotted line) and 
35 Myr (dashed line), all in the center of the initial planetesimal ring. The initial distribution for the simulation is plotted as thin solid line. Note our using 2ns 
instead of * in this and subsequent size distribution plots. Since particles with the size parameter 2ns/ A ~ 1 emit most efficiently, 2ns roughly gives the typical 
wavelength of the emission. This alleviates comparison between the size distribution and the SED (bottom left). Top right: Radial profile of the optical depth 
for the same disk models. Bottom left: Corresponding SEDs. Symbols with error bars are data points (large square, large circle, and large triangle mark 24pm, 
70 pm, and 160 pm excess ratios deduced from Spitzer/MIPS images. Bottom right: Modeled (lines) and observed (symbols) surface brightness profiles at 24, 
70, and 160 pm. The shapes of symbols are the same as in the bottom left panel and in Fig.fT] The shaded areas around the data points indicate the errors. 



the Spitzer point of 0.94 Jy. Let us make a more careful com- 
parison, however. As explained in Sect. 2, the total observed 
flux from 4" outward is unsure, because of an uncertain part 
between 4" and 10". If we consider only the observed flux 
from 10" outward, it reduces to 0.53 Jy. At the same time, 
we can calculate the flux from 10" outward in the reference 
model. The result is 0.40 Jy, which is only slightly below 
the observed one. That value of 0.40 Jy differs insignificantly 
from the calculated flux from 4" outward, 0.43 Jy, because 
the emitting dust in the reference model is entirely located 
outside 10", and it is only the finite PSF's width that "trans- 
fers" 0.03 Jy of the emission closer in. 

As far as the IRAS 25 pm flux is concerned, it is known 
to be quite uncertain due to the large field of view and it is 
inconsistent with the Spitzer 24 pm measurement anyway. 

In the far-IR the model emission matches the Spitzer fluxes 
perfectly and lies within the error bars of the other obser- 
vations. At (sub-)millimeter wavelengths, the measurements 
themselves are contradictory, lying sometimes at more than 
2cr from each other, and it is difficult to judge which of them 
are most accurate. The model SED provides a compromise, 
lying in the middle of the entire set of data points. 

4.5. Radial Surface Brightness Profiles 

Using SUBITO, we calculated the radial surface brightness 
profiles of our disk model and convolved them with the corre- 
sponding PSFs. The final profiles are presented as solid lines 
in the right panel of Fig. [5] The model profiles are not incon- 
sistent with the Spitzer observations. Especially the 70yum 



profile is very close to what was measured. However, both 24 
and 160yum curves are slightly too flat. The 160jt/m profile 
lies under the measurements in the inner part of the disk and 
above them in the outer part, explaining why the total 160 pm 
flux is about right (see the bottom left panel). In contrast, 
most of the 24 pm flux comes from the inner part of the disk 
inside 100 AU. In this region the model profile goes below the 
data points, so that the higher emission in the outer part of the 
disk cannot compensate this deficiency. 

Altogether, we state that the surface brightness profiles are 
more constraining for the disk model than the SED. Already 
our "first-guess" model satisfactorily reproduces the observed 
SED, but the brightness profiles reveal moderate deviations 
from those deduced from the observations. 

5. VARIATION OF MODEL PARAMETERS 

In this section we investigate how the observables (SED, 
brightness profiles) respond to changes in physical parame- 
ters (those of the star, planetesimal belts, and dust, as well as 
the collisional prescription). A specific goal is to check if we 
can improve the agreement of the modeled brightness profiles 
of the Vega disk with observations, preserving the agreement 
in the SED that we achieved in the reference model. Accord- 
ingly, we consider a set of models, the parameters of which 
are listed in Tab. [3] Most of these models differ from the 
reference model by one parameter. We modify several param- 
eters at a time only if these are physically related and this is 
required for consistency. 

For each of the models, we present the results in the same 
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way as for the reference one. The size distribution, optical 
depth profile, the SED, and the radial brightness profiles are 
combined into a single figure (Figs. [6] to [T3T l. each having 
the same structure as Fig. [5] In all the figures, the reference 
model is overplotted with a solid line. The resulting final dust 
masses, final disk masses, and initial disk masses are given in 
Tab. |U 

In the following subsections the variations of the reference 
model are explained and discussed. They are structured ac- 
cording to the underlying physical and astrophysical mecha- 
nisms at work. 

5.1. Delayed Stirring 

Before a debris disk starts to evolve in a steady-state 
regime, a collisional cascade has to ignite and operate for suf- 
ficient time. Ini tiation of the cascade requires a mechanism 
to stir the disk (IWy att 2008). This can be self -stirring by 
largest planetesimals dKenvon & B romlevll2004a|) or stirring 
by planets orbiting in t he inner gap of the disk (IWvatt] [2005b; 
Mustill & Wyatt 2009). External events such as stellar flybys 
can also stir the disk su fficiently, it may hav e been the case 
for Vega ~ 5 Myr ago dMakarov et al.1 120051) . Furthermore, 
after the onset of the cascade, the system needs enough time 
to reach a stead y-state collisional regime at dust sizes (e.g. 
lLohne et al. 2008). Thus the duration of a steady-state disk 
evolution is generally shorter than the system's age. We do 
not know which particular mechanism may have triggered the 
cascade in the Vega disk and how long it is already at work. 
It may have started either shortly after the primordial gas dis- 
persal or much later in the Vega history. 

To investigate the effect of the unknown "collisional age" of 
the Vega disk, we simply took our reference disk model and 
calculated the SEDs and surface brightness profiles at earlier 
time steps. Fig. [5]shows the results at 3.5, 35 and 350 Myr. 
At earlier times, the maximum of the size distribution is 
slightly more pronounced (top left). This results in a moder- 
ate enhancement of thermal emission between 50 and 500 //m, 
which is still in agreement with the observations (bottom left). 
The optical depth profile (top right) shows that 3.5 Myr of col- 
lisional evolution is not sufficient to bring enough particles on 
highly eccentric orbits, so that the profile is steeper than in 
the reference model. The 24 //m and 70 /mi profiles steepen 
(bottom right), the latter being no longer compatible with ob- 
servations. At 35 Myr, the optical depth profile is only shifted 
vertically compared to the reference model, which indicates 
that the spatial distribution has already reached an equilib- 
rium (top right). The final and initial disk masses (Tab.Q are 
close to each other, which is natural as younger disks have 
spent less material in collisions. Besides, the estimated to- 
tal masses of younger disks are smaller than in the reference 
model. The reason is that younger disks that are not in a 
steady-state regi me yet are "dustier " than older disks of the 
same total mass (iKrivov et al.ll2006T) . 

An overall conclusion is that at least several tens of Myr of 
collisional evolution seem to be required to make observables 
consistent with observations. 

5.2. Disk Location 

Our choice of the initial disk extension in the reference 
model comes from resolved images in the sub-mm and radio. 
However, a low resolution of these observations still leaves 
room for reasonable modifications. Hence we now vary the 
initial semimajor axis range of planetesimals, intentionally 



pushing them to the limits posed by the images, in order to 
see the effects more clearly. 

We start with the inner disk edge. By placing additional ma- 
terial closer in, one may expect to increase the warm emission 
and prevent the brightness profile from dropping off towards 
the star too early. Thus, we try shifting the inner edge in to 
50 AU. For completeness, we also add the case with the inner 
edge at 100 AU. The results are presented in Fig. [6] 

As expected, taking ai nne i = 50 AU slightly shifts the SED 
to shorter wavelengths and strongly depresses the sub-mm 
emission. Taking ai nne i = 100 AU yields the opposite, al- 
though the effect is weaker. This traces back to the under- or 
overabundance of larger grains, hundreds of micrometers in 
size (top left). The reason for that, in turn, is the existence of 
two distinct dynamical regimes for bound dust grains. Large 
grains with no or little response to the stellar radiation pres- 
sure essentially inherit their orbital eccentricities from the par- 
ent bodies. For small, barely bound grains, radiation pressure 
is the dominant effect, pushing them to wide, highly eccen- 
tric orbits. Shifting the disk further in increases the average 
collisional velocities, lowers the collisional lifetime of larger 
grains - but not of smaller ones. As a result, the relative abun- 
dance of larger grains is reduced. 

In terms of the surface brightness profiles, it is mostly the 
24 /mi profile that is affected. It rises significantly inside 
200 AU, reproducing perfectly the observations. However, 
in the outer part of the disk it becomes flatter so that the emis- 
sion here is overestimated. The 70 /mi profile remains almost 
unaffected, except for the inner part within 150 AU, which 
responds to fli nne r in the same manner the 24 /mi does, albeit 
less strongly. Finally, the 160 /mi profile preserves its slope, 
but shifts downwards (ai nner = 50 AU) or slightly upwards 

(Oinner =100 AU). 

Similar to the inner edge, an inward shift of the outer edge 
would lower the amount of cold dust, enhancing the warm 
emission. So we changed a out er to 100 and 150 AU to find sim- 
ilar modifications in the dust distribution and thermal emis- 
sion as above (Fig. [7]). Decreasing a outer makes the ring nar- 
rower and shifts the bulk of the material closer in. The max- 
imum in the size distribution becomes more pronounced and 
shifts to smaller grains. The entire SED slightly shifts to- 
wards shorter wavelengths. The peak of the optical depth 
profile becomes stronger and moves closer to the star. This 
directly translates to the radial surface brightness profiles, es- 
pecially at 24 /mi, which becomes appreciably steeper. In- 
creasing fl outer naturally has opposite effects. 

On the whole, it seems that shifting the inner edge of the 
belt inward has a clear potential of getting the 24 /mi profile 
that would better match the observed one. However, the shift 
from 80 AU down to 50 AU that we have tested may be too 
strong, because it may contradict to the sub-mm images. 

5.3. Dynamical Excitation 

We now consider the dynamical excitation of the disk, as 
parameterized by the maximum orbital eccentricities e max that 
planetesimals had at the onset of the collisional cascade. This 
value does not change considerably in the course of the sub- 
sequent evolution (under the assumptions of our model, e.g. 
without planets), and it is approximately the same for disk 
solids of all sizes except for the smallest dust particles that are 
vulnerable to radiation pressure. From the dynamical point of 
view, higher eccentricities increase the collisional velocities 
(altho ugh the collision al rates remain nearly the same, see, 
e.g. lOueck et all2007h and thus the efficiency of the colli- 
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sional cascade. 

We ran two simulations: one with reduced (e max = 0.1) and 
one with increased eccentricity (e max = 0.3). The inclination 
of the disk was taken to fulfill the energy equipartition con- 
dition, / max = e max /2. Since changing e max , but maintaining 
the same distribution of semimajor axes would change the ra- 
dial extension of the disk, we chose to alter aj nner and a out er in 
such a way as to preserve the radial extension of the reference 
model disk. (Strictly speaking, all this applies to the initial 
disk, because we have no control over the disk extension at 
later times in the course of its dynamical evolution.) 

The size distribution in Fig. [8]shows that higher eccentrici- 
ties lead to a depletion of larger grains (> 30 pm) and in return 
to an overabundance of smaller grains close to the blowout 
size. Consequently, the SED drops beyond w 200 /mi and 
rises at shorter wavelengths becoming more narrow. In con- 
trast, with lower eccentricities more grains of radii > 100 /mi 
survive, so that fewer particles with 15 /mi < s < 100 /mi 
are created in collisions, which results in a more pronounced 
maximum between the blowout and 30 /mi. This yields an en- 
hancement of the radio emission and an overall flatter shape 
of the SED. 

The slope of the optical depth profile in the outer part of the 
disk remains nearly unchanged, it just shifts vertically. How- 
ever, in the range of the birth ring the optical depth profiles 
for different e max are different: the lower e max , the broader 
the maximum. It is because for higher e max , the distribution 
of semimajor axes is narrower and the collisional production 
of dust near the center of the ring is much higher than else- 
where. For lower <? max , the semimajor axes are distributed 
more broadly and dust is collisionally produced at compara- 
ble rates everywhere in the ring. 

The brightness profiles respond to the changes in the disk 
excitation in a similar way. In the outer part they just shift 
vertically, and most of the changes are in the birth ring region. 
When e max is reduced to 0. 1 , in the outer disk the 24 //m profile 
matches the observed profile closely. In the region of the birth 
ring it is still too low, as is the reference profile. However, the 
emission now keeps rising inward down to 80 AU, i.e. all the 
way through the birth ring, as does the observed emission. 

5.4. Stellar Luminosity 

Changing the stellar luminosity has a two-fold effect on the 
results. First, it alters the y6-ratio of the dust grains, affect- 
ing their dynamics. Second, a different luminosity affects the 
temperature of the dust grains, thereby changing the SED and 
brightness profiles. Note that all these changes influence only 
the dust portion of the disk, not the larger objects. 

As mentioned above, Vega is a rapid rotator and so the ra- 
diation flux emitted from its surface varies with the stellar 
latitude. In the reference model, we adopted the "equato- 
rial luminosity", 28 L Q . However, the dust disk "sees" not 
only the stellar equator, but also receives stellar radiation from 
higher altitudes. Thu s, we now test the averag e luminosity of 
37 L as derived by Aufdenberg et al. (2006) and, as an ex- 
treme case, the canonical polar value of 57 L (used by many 
modelers before). 

As the luminosity gets higher, the blowout size increases 
(Fig. [3] right), reaching 10/mi for a 57 L central star. The 
entire size distribution shifts horizontally towards larger sizes 
and the jump at the blowout radius becomes more abrupt 
(Fig. [9j. The optical depth profile preserves its shape, but 
moves downward. The reason for the decrease of the optical 
depth level at higher luminosities is simply the increase of the 



grains' blowout size, so that there are fewer grains that could 
stay in bound orbits. 

Interestingly, the SED in Fig. [9] (bottom left) does not 
change substantially, as one would expect from dramatic 
changes in the size distribution. The reason for moderate 
changes seen in the figure can be found by analyzing the size 
distribution (Fig. [4] top right) and the temperature plot (Fig. [4] 
top left). Since for a larger L the size distribution is shifted 
to larger particles, the temperature range of the smallest parti- 
cles, which affects the outer part of the disk the most, becomes 
narrower. Consequently, the SED slightly narrows, too, and 
the maximum at ~ 100 //m becomes more pronounced. 

The changes in the surface brightness profiles are two-fold. 
While the 24 /mi profile steepens with increasing luminosity, 
the 160 and especially the 70 /mi curves flatten. The explana- 
tion for this behavior is as follows. The farther out from the 
star, the faster the temperature decreases with increasing grain 
size (Fig. |4] top left). A comparison with the position of the 
maxima in the size distribution at different distances and for 
different luminosities (Fig. |4] top right) demonstrates that an 
average temperature in the region of the birth ring increases 
with increasing L. In the outer disk the effect is reverse: far 
from the star, a higher stellar luminosity lowers the typical 
dust temperatures. These effects explain why for higher lumi- 
nosities the mid-IR emission rises in the inner disk and drops 
in the outer one, steepening the 24 /mi profile. At the same 
time, the far-IR emission becomes more efficient farther out, 
which flattens the 70 and 160 /mi profiles. 

Our analysis definitely favors an intermediate value of the 
Vega luminosity, exemplified by L — 37 L Q in our tests. 
First, this choice is well justified physically. Indeed, dust 
is exposed to stellar light coming from a range of latitudes, 
thus the "right" luminosity should be between the equatorial 
and polar one. Second, it does provide a better agreement 
with observations. Changes in the 70 and 160 /mi profiles are 
only marginal, so that they still match the observations well 
enough, while the 24 pm profile steepens inside « 250 AU, 
coming much closer to the observed profiles. 

5.5. Chemical Composition 

Like the stellar luminosity, the chemical composition of 
grains affects both the yS-ratio (through the radiation pressure 
efficiency and bulk density) and dust temperatures (through 
the absorption efficiency). 

Mid- and far-I R spectra of some debris disks rev eal dis- 
tinctive features (lJura et al.1 120041: iChen et all 120061) . which 
allows one to get insight into the mineralogy of the dust 
grains. For example, spectra of several disks were matched 
by a mixture of amorph ous and crystalline silicates, silica, 
and s e veral other species dSchutz et al]|2005t iBeichman et al.l 
120051: iLisse ;et al .1120071 12008b . including possibly water ice 
dChen et al. 2008|). Unfortunately, the spectra of the Vega disk 
(available in the Spitzer archive) do not exhibit unambiguous 
features, which poses no observational constraints on its com- 
position. 

In the reference model we used pure astronomical silicate. 
Now, to test possible effects of chemical compositi on, we con- 
sider an astrosilicate matrix w ith water ice (Warren 1983) and 
iron dLvnch & Hunter! [199 11) inclusions. The refractive in- 
dices were calculated according to the Maxwell-Garnett the- 
ory. The amount of inclusions was limited to 10 %, which is 
an upper limit for whic h the effective medium theory sti ll pro- 
vides accurate results (Kolok olova & Gu stafson 200l]). The 
resulting bulk densities for the mixtures with ice and iron 
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FIG. 8. — Same as Fig. [5] but for dynamically less (dashed lines) and more (dotted lines) excited model disks. 
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FIG. 9. — Same as Fig. [5] but assuming a more luminous central star with luminosity of 37 (dashed lines) and 57 L Q (dotted lines). 



The Debris Disk of Vega 



13 



are 3.062 g/cm 3 and 3.757g/cm 3 , respectively (compared to 
3.3 g/cm 3 for pure astrosilicate). The effect of inclusions on 
the absorption efficiencies (except for small grains in the near- 
IR) and /?-ratio (Fig. [3) is minor. Most of the difference in the 
size distributions (Fig. [TUJ top left) probably comes from the 
changes in the bulk density. The blowout for grains with iron 
inclusions is slightly shifted to smaller sizes und the whole 
distribution is stretched, while in the case of water inclusions 
the opposite is true. The radial distribution of dust (top right) 
remains virtually the same. 

In terms of thermal emission, the influence of inhomogene- 
ity is very weak, too. Consistently with the modifications in 
the size distribution, the SED is narrower for ice and slightly 
broader for iron inclusions. The surface brightness profiles re- 
tain their overall shape. Only the 24 fim emission in the outer 
disk becomes slightly stronger and more gently sloping, when 
iron inclusions are present. 

Our conclusion is that inclusions at a 10 % level have only 
minor effect on the observables. We cannot exclude that more 
radical changes in the composition would affect the results 
substantially, but there is currently no observational evidence 
that would justify such changes. 

5.6. Cratering Collisions 

We turn to an analysis of the underlying collisional model 
implemented in ACE. The detailed physics and outcomes 
of binary collisions under the conditions of debris disks are 
poorly known, which represents one of the major sources of 
uncertainty in our simulation results. Accordingly, in this and 
subsequent sections, we vary three key parameters that con- 
trol the treatment of collisions. 

We first explore a hypothetical collisional cascade, in which 
only disruptive collisions operate and the cratering collisions 
do not occur. This means that we only consider collisions 
with specific impact energies above the threshold value (Gjp, 
which shatter both colliders completely. All collisions at 
lower energies (that would in reality erode one or both of the 
colliders) are simply ignored. 

Thus, an efficient way of eroding larger objects by colli- 
sions with much smaller grains is switched off. As a result, 
grains with 20 fim < s < 300 fim are more abundant than in 
the reference model (Fig. ITTb . And conversely, the number 
of grains with s < 20 fim decreases, so that the maximum of 
the size distribution is now effectively shifted to about 30 fim. 
The explanation is simple. Excluding cratering collisions pro- 
longs the collisional lifetime of larger grains, because smaller 
impactors that cannot disrupt but w ould efficiently erode 
them, now leave them intact (see, e.g., Thebault & Augereaul 
2007, and refer ences therein). 

The change in the size distribution shifts the SED towards 
longer wavelengths and slightly narrows it. There is now a 
lack of emission in the mid and far-IR up to lOOyt/m and an 
excess emission at sub-mm. The resulting SED clearly vio- 
lates the data. 

Given the deficiency of small particles on highly eccentric 
orbits, it is not surprising that the optical depth profile be- 
comes very steep. As a consequence, the 70 fim brightness 
profile steepens significantly. In addition, all three profiles 
fall too low. Like the SED, they are no longer consistent with 
the observations. 

We conclude that cratering collisions cannot be ignored. 
They seem mandatory to reproduce the observations of the 
Vega disk with collisional simulations. 



5.7. Energy Threshold for Fragmentation 

In this subsection, we explore the role of the (unknown) 
tensile strength of the solids, parameterized by the shattering 
energy in the strength regime. To this end, we decrease 
Qd, s and b s in Eq. ((TJ in one simulation ("weak material") 
and increase them in another one ("hard material"). 

Figure QT| shows that larger grains benefit from an increase 
of the energy threshold in a similar way they do from neglect- 
ing cratering collisions or from lowering the average impact 
energies. Their collisional lifetime becomes longer and their 
amount increases. And conversely, a lower reduces the 
amount of larger particles. For smaller grains, the number 
of potentially hazardous impactors is not determined by the 
change in the critical impactor mass for disruption (that comes 
along with the change in critical energy) but by the blowout 
limit, which remains unchanged. Consequently, decreasing 
(2p "supports" smaller grains, producing a more pronounced 
first maximum in the size distribution. 

In the SED, a lower critical energy leads to a strong shift 
of the maximum to about 80 //m, and makes the rise in the 
mid- to far-IR steeper, whereas the sub-mm and millimeter 
part lowers and flattens. The opposite changes, albeit less pro- 
nounced, are seen for a higher critical energy. In both cases, 
the agreement between the modeled and observed SED be- 
comes rather worse. 

Similarly to the size distribution, the optical depth profile 
responds to a harder material in nearly the same way as to ex- 
cluding the cratering collisions. As far as the surface bright- 
ness profiles are concerned, the only real improvement can 
be found in the 160pm profile for larger Qd s and b s . How- 
ever, this is accompanied with a steepening and flattening of 
the 70 and 24 fim profiles, respectively, which are then clearly 
inconsistent with the observations. 

We conclude that using "weak" or "hard" material with re- 
spect to the nominal one does not generally improve agree- 
ment with the observations. When improving one of the three 
surface brightness profiles, for instance, this makes one or two 
of the others worse. We find, however, that results are very 
sensitive to the critical energy. Thus moderate modifications 
in the critical energy can be useful for "fine-tuning" of the 
models. 

5.8. Fragment Distribution 

One more essential part of the collisional description is the 
distribution of fragments produced in a single collision. In the 
reference model we assumed their mass distribution to follow 
a pow er-law with an index 77 = 1.833 retrieved from experi- 
ments (Fui iwara et al.|[l977l) . but the experimental conditions 
do not necessarily reproduce the conditions of debris disks. 
Here, we try another two disk models, one with an enhanced 
production of small particles (77 = 1.95) and one with a re- 
duced production (77 = 1.6). 

The effect on the size distribution in Fig. [12] is not very 
strong. An increase of 77 enhances the production of small 
particles so that the total distribution becomes flatter. This 
makes the SED broader: the far-IR emission decreases while 
the sub-mm fluxes are enhanced. Reducing 77, however, trims 
the production rate of small particles, so that the maximum in 
the size distribution becomes broader and is shifted to about 
15 fim. Consequently, the SED becomes somewhat narrower, 
with a steeper rise in the mid-IR, stronger emission in the far- 
IR, and a steeper fall-off in the sub-mm. These changes are 
minor, so that the SEDs for both 77 values are consistent with 
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the observed SED, as is the SED in the reference model. 

Changes in the optical depth are subtle, but not unimpor- 
tant, particularly at the outer edge of the birth ring. While a 
lower j] makes the optical depth profile smoother, a higher 
77 creates a slight dip at the outer end of the planetesimal 
belt. In this radial zone, emission stems predominantly from 
intermediate-sized grains, which are placed by radiation pres- 
sure in moderately eccentric orbits, but still cannot reach 
the outermost regions of the disk. Reducing 77 increases the 
amount of these particles compared to the reference model, 
which smoothens the optical depth profile in this region. 
And conversely, an increase of 77 depresses the population 
intermediate-sized grains, causing the dip. 

In the radial surface brightness profiles these changes are 
evident in the 24 pm profile, which becomes flatter for small 
77 and steeper for high 77 in the region up to about 200 AU. The 
vertical shifts are in agreement with the modifications in the 
SED. 

A conclusion is that a flatter size/mass distribution of frag- 
ments with 77 = 1 .95 makes the 24 pm brightness profile more 
consistent with the observed ones, without making other two 
profiles and the SED worse. 

5.9. The Best Fit 

Different modifications in the previous section have shown 
no simple way of further improving the agreement of our ref- 
erence model with the observations. However, we found that 
variation of some parameters is able to change the results in 
the desired direction. We now combine several of the modifi- 
cations that looked promising: the disk was extended inwards, 
the eccentricities were reduced, the luminosity increased, and 
a steeper fragment distribution was assumed. Specific param- 
eter values are listed in the last line of Tab. [3] 

The result is depicted in Fig. [T3]with dashed lines; as al- 
ways, solid line shows the reference model for comparison. 
In terms of the SED (bottom left), all photometry data (aside 
from the IRAS 25 pm points) are reproduced within the error 
bars shortward of 200 pm. At longer wavelengths, the ob- 
servational data themselves split into groups that are not in 
agreement with each other. The model perfectly matches the 
upper set of points. 

In terms of the radial surface brightness profiles, the 160 /mi 
profile is nearly the same as in the reference model. The 
24 /mi profile is at about 1 <x in all regions of the disk. From 
all three curves the 70 /mi profile is the closest to the observed 
one. 

6. DISCUSSION 

6.1. Modeling the Vega Disk "from the Sources" 

In this paper we have attempted to reproduce observations 
of the Vega disk by assuming that the observed emission 
stems from the debris dust, which is produced by a "Kuiper 
belt" at ~ 100 AU in a steady-state collisional cascade. To 
this end, we have performed involved collisional simulations 
to generate and evolve the disk of solid material from plan- 
etesimals to dust and then calculated thermal emission of the 
dust portion of the disk to confront the results with the avail- 
able observations. Compared to the commonly used modeling 
technique, in which dust distributions that provide the best fit 
to the data are sought, we seek parent planetesimals whose 
collisional debris could produce emission that matches the ob- 
servations. This is obviously one major step more than in the 
traditional method: we try to go back to sources by modeling 



actual physical processes that operate in a debris disk. 

Starting this research, we wanted to test whether our ap- 
proach can withstand a detailed comparison with various sets 
of observational data on one particular, and well studied, de- 
bris disk system. We have chosen Vega, an archetype debris 
disk. One particular motivation for this choice was that previ- 
ous studies uncovered peculiarities, casting doubt on a "stan- 
dard collisional cascade" as the main mechanism sustaining 
the disk of Vega. 

We understood from the very beginning that our efforts may 
not be rewarding. The fact that our approach involves a phys- 
ical modeling to the whole "depth" automatically adds major 
uncertainties, notably those related to the collisional physics. 
Contrary to the standard data fitting, we cannot control size 
and radial distribution of visible dust and we are not at lib- 
erty to add, for example, another population of dust to im- 
prove agreement with observations. Thus it would be unreal- 
istic to expect that our approach could immediately deliver re- 
sults that are consistent with observations to that same extent 
as in the traditional method. It was to our surprise therefore 
that already with a 'first guess' reference model we are able 
to reproduce the SED and the modeled mid- to far-IR radial 
brightness profiles were not inconsistent with Spitzer/MIPS 
data. We were then able to further improve the agreement by 
variation of model parameters. 

6.2. Blowout or bound grains? 

It is interesting to trace, why lSu et alJ (120051) needed an ex- 
tremely high amount of blowout grains to explain the MIPS 
observations, whereas we are able to reproduce the same data 
without any blowout grains (in our models, they make a neg- 
ligible contribution to the SED and radial brightness profiles 
at all wavelengths considered and at any distance in the disk). 
By fitting the MIPS photometr y and radial pro files with a sin- 
gle power-law size distribution. lSu et 51(120051) found the best 
fit to be s m i n = 1 /mi and s max = 50 pm, with a slope of -3.0. 
Assuming nominal stellar luminosity of L — 60 L Q , already 
the compact grains have a blowout limit of a: 8 pm. How- 
ever, if the grains are highly porous, which of course they 
may, then the blowout limit will shift to much larger sizes, 
so that all grains in the 1 /mi... 50 /mi could indeed be in un- 
bound orbits. An additional argument to favor the radiation 
pressure-induced outflow seemed to be the deduced bright- 
ness profile at 24 pm with a slope of -3... -4, since a -3 slope 
is what blowout particles with nearly constant "terminal" ve- 
locities would produce. At that time, it was not yet realized 
that slopes in the range -3. ..-4 would equally be typical of an 
extended disk of small bound grains in elliptic orbits around 
a parent planetesimal ring, as was found later numerically 
dKrivoy et alj|2006h It hebault & Augereaull2007l) and analyt- 
ically dStrubbe & Chiangll2006l) . 

In this study, we assume compact grains and L — 28 L 
in the reference model (note that the very idea to reduce 
the stellar lumi nosity would sound strange at the time of the 
(Su et al. 2005) study, because the fast rotation of Vega was 
not yet discovered). A s a resu lt, the blowout radius reduces 
to as 4-pm. Thus where ISu et alJ had 1 /mi... 50 /mi grains, we 
have 4pm.. .50pm ones. For a -3.0 size distribution slope, 
the emitting cross section area is equally distributed over 
the sizes. Thus, even if we take into account that smaller 
grains are somewhat hotter than larger ones, by excluding the 
1 /mi. ..4 pm subra nge we d o not lose much of the 24 pm emis- 
sion compared to ISu et alJ The radial distribution of dust in 
our model, as explained above, is not very different either. 
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FIG. 12. — Same as Fig. [3] but for model disks with a flatter (dashed) and steeper distribution of fragments (dotted). 
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That is why we arrive at a similar level of 24 /zm emission. 
In terms of dust mass, the difference is even smaller. Indeed, 
the dust masses we derive here (6.6 x 10~ 3 M e in t he refer- 
ence model, see Tab.HI are close to those derived by lSu et all 
(2.8 x l(T 3 M e ). 

6.3. Mass loss from the disk 

We now make estimates of the mass loss from the disk. 
The total mass of dust (up to 1mm) in the reference model 
is ~ 7 x KT 3 M e (TableO. Assuming for simplicity that the 
size distribution follows a -3.5 power law, the mass of the 
smallest bound grains (say, up to 10/zm) is ~ 7 x 10~ 3 M© x 
-y/lOjum/lmm ~ 7 x 10~ 4 M e . The steady-state mass of 
blowout grains is then by a factor of 100 smaller (a strength 
of the dip in the size distribution, which is the ratio of the 
collisional lifetime of bound grains to the disk-crossing time 
of unbound ones, see Fig. |4] right), giving ~ 7 x 10~ 6 M ffi . 
Their lifetime is ~ lOOOyr, so that the mass loss rate is 
~ 7xlO" 9 M e yr 1 . Therefore, over the system's age, 350 Myr, 
the disk must have lost ~ 2M e of material. This estimate is 
consistent with the difference between the initial and final disk 
masses given in Tab. 2] typically a few M @ . 

6.4. The 24 pm Emission 

It was "the 24 /zm problem" - an apparently too strong and 
radially extended 24 /zm emission compared to what was ex- 
pected from dust in bound orbits — that triggered debate on 
whether the Vega disk contains an excessive number of small 
blowout grains, in compatible with a steady-state collisional 
cascade dSu et al J 120051) . Thus we now discuss in more de- 
tail how the 24 /zm flux predicted by our models compares to 
Spitzer data. 

Although the models presented here are in a reasonable 
agreement with observations, most of them somewhat under- 
estimate the observed 24 /mi emission in the parent ring re- 
gion, at 80-120 AU or 10"— 15", while slightly overestimating 
it farther out from the star. As a result, the total flux outside 
10", which is dominated by the flux from the ring region, is 
slightly below the observed value. For example, our reference 
model predicts 0.40 Jy outside 10", while the observed flux is 
0.53 Jy. These deviations are subtle and probably not of seri- 
ous concern. We argue that they may s imply be caused by the 
roughness of the model (see Sect 16.71 . Indeed, w e we re able 
to find a combination of model parameters (Sect. 15.91 ) which 
reproduced the observed 24 /mi profile outside 10" quite well. 

However, the question of 24 /mi emission from the inner 
system (< 10") remains open. Our analysis in Sect. |2] yields 
a total 24 /zm flux from 4" (30 AU) outward of 0.94 Jy, al- 
though the true value may be lower, because central part of 
the MIPS images is saturated. Assuming, however, the data 
to be accurate, the rise of the 24 /zm flux from 10" inward can 
hardly be explained with the models presented here. A natural 
explanation would be an additional dusty belt in the system at 
~ 10 AU. Such a belt could enhance the 24 /an emission com- 
ing from the "main" disk. As yet it is not clear, however, if 
any constraints on such a belt can be found in the Spitzer/IRS 
spectrum of Vega. Nor is it clear whether the inner system 
may accommodate such a belt if, as conjectured, it hosts one 
or more close-in planets. In the future, this simple hypoth- 
esis could be checked or falsified directly, for instance with 
mid-IR interferometry. 

6.5. The 850 /mi Emission 



So far, we confined our analysis to resolved images in the 
IR and did not consider explicitly sub-mm and radio images. 
The main reason for that is a low resolution of such measure- 
ments, implying that only weak constraints can be put on the 
radial brightness profiles at long wavelengths. However, at 
least the total flux at sub-mm wavelengths derived from the 
images serves as an additiona l test t o the models. 

In an analysis by ISu et al. (2005), the observed sub-mm 
emission could not be reproduced with a two-component dust 
disk (2 /mi and 18 /mi in radius) that was sufficient to fit all 
available data at shorter wavelengths. To cope with the prob- 
lem, they artificially added a population of larger grains, with 
a radius of 215 /zm. In our approach, solids from dust to plan- 
etesimals have a continuous size distribution, which is not 
postulated, but physically modeled. From Figs. [514X31 (left 
bottom panels) it is apparent that our simulations naturally re- 
produce the sub-mm flux with a reasonable accuracy. An ad- 
ditional consistency check that we made was to calculate the 
850 /zm profile of our best fit model (Sect. IBT91 and convolve it 
with a Gaussian of 16" beam size. The resulting pro file was 
then compared with the SCUBA profile extracted by Su et alJ 
(2005) from the original images published by (Holland et al.1 
1998). The modeled profile is slightly narrower than that ob- 
served and the maximum between 50 and 100 AU is by a fac- 
tor of two lower. Given the large width of the PSF and cal- 
ibration uncertainties of SCUBA observations, and that Mie 
calculations lik ely u nderestimate sub-mm emission (as dis- 
cussed in Sect. 16.71 . we deem the agreement with the data 
satisfactory. 

6.6. Model Parameters 

As we saw, even our "best fit" (Fig. [T3l ) cannot match the 
observations perfectly. The discrepancies could stem either 
from the fact that our "first-guess" choice of parameters in 
the reference model was not the best, from the limitations or 
shortcomings of the collisional and thermal emission model, 
or for the two reasons together. In this and the next section, 
we address the both possible reasons in turn. 

In Sect. [5] we investigated in detail how the dust distribu- 
tions and thermal emission are affected by a large array of 
physical parameters: 

• The "collisional" age of the disk, i.e. the time elapsed 
from the onset of the steady-state collisional evolution 

• The location and extension of the parent planetesimal 
belt 

• The dynamical excitation of the belt, parameterized by 
the maximum orbital eccentricity of planetesimals 

• The stellar luminosity 

• The chemical composition of the visible dust 

• Several parameters that control individual collisions be- 
tween the objects in the disk, from planetesimal to dust. 
In particular, we checked the role of cratering colli- 
sions, the critical impact energy threshold for disrup- 
tive collisions, and the distribution of fragments after a 
collision. 

We have identified two parameters that have a major influ- 
ence. One is the stellar luminosity, and its effect is particularly 
interesting. One might think that a more luminous central star 
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would increase dust temperatures, but this is untrue. A higher 
luminosity implies a larger blowout size. Thus the most abun- 
dant grains in the disk — those just above the blowout limit — 
are now larger and therefore colder (Fig.|U top left). The net 
result is that, counterintuitively, the characteristic dust tem- 
perature in the disk does not change much (in the outer disk it 
even decreases) with increasing luminosity of the central star. 
Another parameter is the efficiency of cratering collisions. In 
the extreme case where these are switched off, the simulation 
results contradict the observations. 

Contrary to our expectations, it turned out that other param- 
eters probed have only minor to moderate effect on the SED 
and surface brightness distribution of the disk, with conse- 
quences being two-fold. On the one hand, this implies that 
the model predictions are rather robust. As an example, our 
reference model was already in a rough agreement with ob- 
servations, so that further improvements through "parameter 
tuning" seemed easy, but these turned out to be difficult. We 
find meaningful parameter sets that provide better fit to obser- 
vations than the reference model does (see Fig. [13}, but the 
agreem ent is somewhat worse than in the study of ISu et aT] 
(2005). On the other hand, a weak dependence of the ob- 
servables on many parameters restricts the possibility of con- 
straining them, which is somewhat unlucky, because placing 
constraints on model parameters is one of the important goals 
of the modeling. 

6.7. Limitations of the Model 

We are fully aware that our modeling approach, as every 
other, involves a number of simplifying assumptions that may 
limit its applicability and influence the results. Here, we list 
the most important caveats. 

Many assumptions have been made in describing collisional 
physics. Our collisional prescription approximates the criti- 
cal shattering energy with two power laws (Eq.[T}, which may 
be pa rticularly crude at dust sizes (e.g. Thebault & Augereau 
2007). The mass of the largest fragment and the distribu- 
tion of smaller debris may deviate from what was assumed 
here, and any real disk should be composed of objects whose 
mechanical properties (and even the bulk density) vary from 
one object to another (e.g. pre-shattered objects could be less 
dense and more loosely bound than pristine ones). 

A major simplifying assumption in treating the dynamics of 
planetesimals and their dust is that we ignore alleged plane- 
tary perturbers int erior to the main belt. The consequences are 
discussed in Sect. 16.91 At dust sizes, we di d no t take into ac- 
count P-R drag, its role is discussed in Sect. 16. 81 In calculating 
the radiation pressure force acting on dust grains, we assumed 
them to be compa ct and spherical, th us ignoring possible non- 
radial effects (e.g. lKimura et al.l2002l) . Like mechanical prop- 
erties, optical properties of dust may vary from one grain to 
another, resulting in different response to radiation pressure, 
different temperatur es, and different th ermal fluxes even for 
like-sized particles (iKrivov et al.l 12006). Furthermore, even 
spherical particles are treated in an approximate way. We 
applied Mie theory to model the emission properties. Al- 
though this method is classical and commonly used, it should 
be treated with caution. One particular concern is that Mie 
calculations probably underestimate the emi ssion in the sub- 
mm and radio due to neg lected size effects ( Stognienk o et alj 
[1993 IKrivov etd]l2008D . 

6.8. The Role of the Poynting-Robertson Effect 



Poynting-Roberston (P-R) drag mostly affects smallest par- 
ticles and thus emission at shortest wavelengths considered. 
The P-R force moves such grains inward, placing some of 
them interior to the inner edge of the birth ring. The warm 
emission of these particles especially around the inner edge 
of the birth ring should increase. 

However, with the Vega disk's relatively high optical depth 
(8.3 x 10~ 4 at 100 AU in the reference model), the colli- 
sional timescales of dust grains are shorter than timescales 
over which Poyining-Roberston (P-R) drag causes their ap- 
preciable radial displacement (Fig. [2}. Thus the Vega disk can 
be referred to as a colhsion-dominated, rather the n transport- 
dominated disk (IKrivov et al.ll2000t IWvattll2005al) . Still, it is 
useful to check to what extent PR drag may affect the results 
in terms of SED and brightness profiles. 

To this end, we have switched on P-R drag in our reference 
model. Unfortunately, when a drag force, the P-R force in our 
case, is a dded to a collisional model, the mass-time scaling 
law (Sect. 13.11 ) is no longer valid. One has to take the initial 
disk mass that yields a disk with the "correct" dust mass af- 
ter 350 Myr, i.e. the dust mass that gives the maximum of 
the SED at the level actually observed. This means trial and 
error, i.e. several ACE runs with different initial disk masses 
followed by SEDUCE and SUBITO runs. What makes the 
modeling even more demanding, is that the presence of a drag 
force implies diffusion in the phase space of pericenters and 
eccentricities, which slows down each ACE run appreciably. 
We had to perform four ACE runs, each of which took about 
20 core-days CPU time. 

The "right" dust mass after 350 Myr of evolution with P-R 
was achieved when the initial disk mass was set to 20.5 M e 
(instead of 18.9 M e in the reference model without P-R), and 
the final disk mass was 18.0 M e (instead of 16.3 M e without 
P-R). As expected, the influence of P-R on the 160 /im and 
70 //m turned out to be completely negligible. At 24 /mi, the 
emission in the outer disk increases by ~ 10 % and in the birth 
ring (at 80 AU) by ~ 60 %. Thus the whole 24 fim profile gets 
somewhat steeper, and agrees with observations slightly better 
than the original profile in the reference model (solid lines in 
Figs. I5UT31. However, the improvement is only minor, and 
we conclude that the P-R effect can safely be neglected in 
modeling the Vega system. 

6.9. Presumed Planets in the Vega System 

One major caveat not discussed in Sect. l6.7l is that our colli- 
sional model, implemented in the ACE code, ignores effects of 
a possible planet (or planets) interior to the planetesimal belt. 
Below we briefly outline the facts that point to the presence 
of such planets in the Vega system and discuss how these per- 
turbers may affect the observed properties of the debris disk. 

Asymmetries in t he Vega disk were first discovered by 
iHolland et ail (119981) in a SCUBA 850//m image, and sub- 
sequent sub-mm and ra dio observat i ons ha ve confirmed a 
clumpy ring structure. Wiln er et al.l {2002) introduced the 
idea of a Jupiter mass planet trapping dust in mean-motion 
resonances. They applied A^-body simulations and thermal 
emission calculations to model this scenario and achieved a 
reasonable agreement with their IRAM map. An in-depth 
in vestigat i on on the Vega system dynamics was performed 
by I Wyattl Q2003) who suggested that a Neptune-mass planet, 
migrating outward from 40 to 65 AU over a time span of 
~ 56 Myr, may have cleared the inner part of the assumed 
planetesimal disk and trapped a significant amount of mate- 
rial in the 3 : 2 and 2 : 1 resonances, thus creating two 
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clump s as seen bv lHolland et al.l dl998l) . Later on. lReche et alj 
(2008) generalized this theory to account for ec centric plane- 
tary orbits. Their findings are similar to those of lWvattl (2003) 
with the difference that they require a Saturn-mass planet on 
a low eccentricity orbit to account for the brightness asym- 
metry. For the clumps to be visible against the non-resonant 
background, low planetesimal eccentricities of < 0. 1 are nec- 
essary. Planets with masses greater than w 2Mj up i ter would 
raise planetesimal eccentricities to about 0.2. Besides, too 
massive planets would quic kly deplete the disk. 

As shown in section l53l our simulations slightly favor low 
planetesimal eccentricitie s up to 0.1 . This is in agreement 
with the limit given by iReche et all d2008h . Still, a ques- 
tion arises whether non-inclusion of the azimuthal structure 
in our simulations (ACE treats rot ationally-sym metric disks) 
is a reasonable assumption. Indeed, Wyatt (2006) investigated 
the dust production in a clumpy disk of resonant planetesi- 
mals and showed that local dust production from the clumps 
is strongly enhanced and conversely, it is depressed between 
the clumps. We argue, however, that the net effect on the 
SED and radial profiles of brightness is much weaker, be- 
cause these depend on th e azimuthally-averaged dust produc- 
tion rates. lOueck et al.l (120071) found that the average col- 
lisional rate in a resonant planetesimal belt is typically not 
more than twice as high as in a similar non-resonant belt, 
while the average collisional velocities are nearly unaffected 
by the resonant clumping. Thus a collisional cascade in a 
resonant, clumpy belt can be approximated by a cascade in a 
non-resonant, rotationally-symmetric belt with the same mass 
at the same location, but somewhat higher orbital eccentric- 
ities of planetesimals to mimic moderately enhanced colli- 
sional rates. We stress that this is only valid when consid- 
ering azimuthally-averaged observables, not the azimuthally- 
resolved structure seen in the images. For example, our ap- 
proach is not suitable to make predictio ns for spiral s tructure 
expected to emanate from the clumps. IWyattl y006) argues 
that such structure should be seen in mid- to far-IR images, 
and that it is not, may simply be due to insufficient resolution 
of the Spitzer/MIPS images or confusion in the photospheric 
subtraction. 

Throughout the study, we assumed that the initial eccentric- 
ities and inclinations of the parent bodies in the planetesimal 
belt are distributed according to energy equipartition. This 
assumption would be reasonable if the distribution of orbits 
was controlled by mutual collisions and gravitational scatter- 
ing among planetesimals, but it may not hold as soon as reso- 
nant interaction with planets occurs. A well-known example 
is our Edgeworth-Kuiper belt, in which the eccentriciti es and 
inclin ations of objects are distributed differently (e.g. Brown 
l200l . 

Apart from the suspected planet that sculpts the main belt, 
the Vega system may contain more planets closer in. In 
fact, a damped outward migration of the presumed planet 
that explains the clumps requires the presenc e of another, 
more massive planet in the system closer in (e.g. Go mes et alj 
l200l). An inner plan et, or planets, could stir the disk 
(Musti ll & Wvattl 120091) . Furthermore, several planets to- 
gether could produce intricate combined dynamical effects on 
the main planetesimal belt and its dust. However, it seems 
premature to discuss them until new observations have deliv- 
ered evidence for these planets. 

6.10. The Exozodi in the Vega System 



As mentioned in the introduction, dust was surprisingly 
discov ered in the innerm ost part of the Vega system, inside 
1 AU ( Absil et al. 2006). Although reminiscent of the zodia- 
cal cloud of the solar system, this "exozodi" of Vega remains 
a mystery. It seems to be far too dusty, and the grain sizes 
retrieved from observations far too small, to be explained by 
collisions in an "asteroid belt" or evaporation of comets. One 
possibility would be a transport of planetesimals from the 
"main" debris disk inward and their subsequent disruption or 
evaporation. Such a transport would require the presence of 
at least two planets. In fact, a two-planet configuration - a 
"Jupiter" inside and a "Saturn" outside that shapes the main 
disk - could suffice (Vandeportal et al., in prep.). Thus, the 
very existence of the exozodi may strengthen the expectation 
that Vega hosts several planets, as discussed above. 

The direct contribution of the exozodi to the emission at 
24 fim amounts to * 0.6 Jy. This is about the emission which 
is lacking provided the used photometry data is accurate. 
However, as the exozodi could not be resolved with Spitzer 
and the image is saturated at the stellar position, this very in- 
ner part of the Vega disk cannot have affected the observation 
in the outer parts of the system and can therefore be consid- 
ered negligible. Still, if not directly, the Vega exozodi could 
have an indirect impact on the measured dust emission. Dust 
inside 1 AU could have a shielding effect on dust located far- 
ther out. However, a simple estimate shows that the amount 
of stellar radiation to which outer dust is exposed would only 
reduce by a factor of ~ 10~ 5 . Thus, we conclude that the 
very inner part of the system has no impact on the outer disk's 
emission analyzed in this paper. 

7. CONCLUSIONS 

Our analysis suggests that the debris disk of Vega is com- 
patible with a standard scenario, in which visible dust orig- 
inates from a steady-state collisional cascade operating in a 
"Kuiper belt", whose existence at ~ 100 AU from the star 
is evident in sub-mm images. We model the dust production 
from the sources, and find that thermal emission of the result- 
ing dust is fully consistent with the photometric data across 
the entire wavelength range from mid-IR to radio covered by 
observations. Furthermore, we are able to naturally reproduce 
the radial brightness profiles at 24, 70, and 160//m derived 
from Spitzer/MIPS observations. Finally, the modeled emis- 
sion agrees with the low -resolution images at 850 pm taken 
with SCUBA. 

If the Vega disk is maintained by a steady-state collisional 
cascade, which appears likely, its total mass (in < 100 km- 
sized bodies) must fall in the range from several to several 
tens of Earth masses. Provided that collisional cascade has 
been operating over much of the Vega age of ~ 350 Myr, the 
disk must have lost a few Earth masses of solids during that 
time. Further constraints of the parameters of the system and 
physical processes operating in the disk are as follows. A 
reasonable amount of stirring should be present in the plan- 
etesimal ring. We demonstrate that planetesimals are likely to 
have eccentricities of the order of ss 0.1. ..0.3, but the origin 
of stirring cannot be constrained. It may come, for instance, 
from a presumed giant planet interior to the belt, or the disk 
can be self-stirred by largest, Pluto-sized planetesimals. Next, 
we show that the modeling results sensitively depend on the 
luminosity of the central star. Importantly, in the particular 
case of Vega, using a reduced radiation flux from the stellar 
surface at low latitudes, which was derived from its fast ro- 
tation, is mandatory to match modeled dust emission to the 
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data. Another important prerequisite for this is to include cra- 
tering collisions between the disk solids as a collisional out- 
come into the model. 

Another goal of this study was to use the Vega disk as a 
stringent test for our modeling approach, which is a physi- 
cal m odeling of a debris disk "from sources" (iKrivov et al.1 
2008). We deem the test successful. We have shown that the 
approach does work and has a potential to deliver constraints, 
most notably on the properties of directly invisible planetes- 
imals, that are not possible to put with other methods. Our 
collisional and thermal emission model needs to be further 
tested and "calibrated" on other resolved debris disk systems. 
Then, it can be used as a routine procedure in dynamical mod- 
eling of debris disks, both known and expected to be detected 
with facilities like Herschel or ALMA. 

The data reduction and modeling presented here describe 
and explain the outer system, outside m 80 AU. Inside 
that distance, our analysis of 24 //m emission observed by 
Spitzer/MIPS indicates a possible rise of the 24 /mi flux from 
10" inward, although it is not certain, because central part of 
the MIPS 24 fim images is saturated. Assuming, however, the 
data to be accurate, they cannot be explained by the models 
presented here. A natural explanation would be an additional 
dusty belt in the system at ~ 10 AU. Such a belt could en- 
hance the 24 fim emission coming from the "main" disk. In 
the future, this hypothesis could be checked or falsified di- 
rectly, for instance with a more accurate mid-IR photometry 
of mid-IR interferometry. On any account, further observa- 
tional and theoretical effort is required to shed more light onto 
the inner part of the Vega system, including alleged planets, 
possible planetesimal belts and dust rings, of which one — an 
"exozodi" at just ~ 1 AU — has recently been discovered with 
near-IR interferometry. It is the inner system that must even- 
tually bring clues to the entire architecture and the formation 
history of the Vega system. A better knowledge of the inner 
part of the system, most notably suspected planets there, will 
also result in a better understanding of how the outer debris 
disk operates. 
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